Many-body Diagrammatic Expansion for the Exchange-Correlation Kernel in Time 

Dependent Density Functional Theory 
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A diagrammatic expansion for the dynamic exchange-correlation kernel f xc of time dependent 
density functional theory is formulated. It is shown that f xc has no singularities at Kohn-Sham 
transition energies in every order of the perturbation theory. However, it may diverge with the 
system size in extended systems. This signifies that any approximate perturbative substitute for f xc 
requires a consistent perturbative treatment of the response function to avoid uncontrollable errors 
in the many-body corrections to excitations energies. 
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The central problem of time dependent density func- 
tional theory (TDDFT) 1 is to find an adequate approx- 
imation for the dynamic exchange-correlation (xc) po- 
tential v xc - In contrast to static DFT, where the lo- 
cal density approximation (LDA) has been extremely 
successful, no universal recipe for a dynamic v xc has 
been found, and it remains unclear if such a recipe ex- 
ists. The adiabatic LDA (ALDA), which is most popu- 
lar in practical TDDFT calculations, is valid, by its na- 
ture, only in a quasistatic limit. This may suffice for 
a real-time dynamics of melting 2 or desorption 3 , but is 
hardly relevant for electronic excitations in insulators, 
where ALDA predicts the same erroneous band gaps as 
the static LDA. 4 The improvement in excitation spec- 
tra of atoms and molecules which has been obtained in 
TDDFT using ALDA or the optimized effective poten- 
tial (OEP) approximation is not indicative because in 
small systems the correction to the Kohn-Sham excita- 
tion energies is very small and approximations for f xc 
play a secondary role in comparison to the accuracy of 
the static xc potential. 5,6 

In general, electron-hole (e-h) excitation energies can 
be found as the poles of the density response function 
x(r,r',w). TDDFT relates the exact x(w) to the sus- 
ceptibility of non-interacting Kohn-Sham (KS) particles 
Xs(w) via the equation 



(1) 



with V(oj) = Vc + fxc(u), where Vc is the Coulomb 
potential and the xc kernel f xc = Sv xc (r,t)/Sn(r' ,t') de- 
scribes xc effects at the linear response level. 7 Equation 
(1), which simultaneously accounts for the self-energy 
and the electron-hole correlations, may seem like an at- 
tractive alternative to the very laborious two-step ap- 
proach involving a GW calculation for one-particle states 
with a subsequent solution of the Bethe-Salpeter (BS) 
equation for an e-h pair. Unfortunately, no approxima- 
tion for f xc is available which would be as efficient as LDA 
in the static DFT. In ALDA a time-dependent density is 
simply inserted in the LDA xc potential, and f xc becomes 
an instantaneous point interaction, which is qualitatively 
different from a nonlocal and retarded xc kernel in sys- 
tems with an energy gap. In fact, it is a priori not clear 



whether f xc allows any reasonable approximation, since 
its analytical properties in a non-homogeneous system 
are not known. 

From the viewpoint of the standard many-body for- 
malism V(uj) in Eq. (1) acts as a mass operator for 
the density propagator x, similar to the self-energy S = 
Gq 1 — G — 1 , with G and Go being interacting and non- 
interacting one-particle Green's functions. In contrast to 
the Green's function, E has no free-particle singularities 
in any order of the perturbation theory, as it contains 
only irreducible diagrams. Hence it allows perturbative 
approximations. Similarly, the introduction of the "den- 
sity mass operator" V = x^ 1 ~ X 1 is motivated only 
if the latter does not possess the e-h singularities, which 
are present in xs and any finite order approximation to 
X- 

In this paper we use the KS-based many-body dia- 
grammatic technique of Ref. 8 to derive a perturbative 
expansion for f xc . We find that the kernel is indeed reg- 
ular at KS frequencies. Yet, caution must be exercised 
in applying perturbative approximations [such as OEP 
(Refs. 7,9)] for f xc in large systems to avoid uncontrol- 
lable errors or even unphysical divergences. 

We start with the equation which relates x(w) to the 
proper polarizability x 

= xH + xHVbx(w)- ( 2 ) 

It is convenient to split x as X — Xs + n xc where the xc 
part ir xc can be represented as a series of graphs. 8 The 

first-order contribution -k x x } and examples of the Second- 
ly") 

order corrections tt xc are shown in Fig. 1, where solid 
and dashed lines stand for the KS Green's functions and 
the Coulomb interaction respectively. To the wiggled line 
is assigned the inverse KS susceptibility Xg 1 - It describes 
the scattering of KS particles by the xc potential v xc . By 
the definition of v xc , these processes exactly compensate 
the change of the density due to the self energy insertions 
in every order of the perturbation theory. 8 

The graphical expansion of f xc can be obtained from 
the relation f xc = x^ 1 — x -1 (Ref. 7), which follows from 
Eqs. (1) and (2). The formal expansion of this equation 



%:: 

(b) 

— <C_ 



o 



►O-O 



> 



;r->~<)III...n>~ 




FIG. 1: (a) The first-order correction 7ri c to the proper po- 
larizability; (b) examples of the second- order corrections 



in terms of ir xc gives the series 

fxc = X^^xcXs 1 ~ X~s X ^xcX~s X ^xcX~s X ■■■ (3) 

Inserting ir xc — n^J + t^xI + t^xc + • • • and collecting 
the diagrams of the same order we obtain fx]) , /ic\ etc. 
The n-th order contribution f$ is a sum of the polariz- 
ability of the n-th order irxc with two attached wiggled 
lines and all diagrams with lower-order polarization loops 
connected by wiggled lines (Fig. 2). A closer inspection 
shows that all graphs with internal wiggled lines can be 
generated from a limited set of graphs with no internal 
but two external wiggled lines. This leads to the fol- 
lowing rules for constructing f^c'- (i) Draw all standard 
diagrams 10 for the proper polarization operator of the 
n-th order and attach wiggled lines to external points 
(construction of parent graphs), (ii) If possible, separate 
any given graph into two by cutting two fcrmionic lines. 
Join the external fermionic lines of these parts, connect 
them by a wiggled line and change the sign. Do not cut 
lines attached to the same wiggled line, (iii) Apply (ii) 
to all possible cuttings in all graphs, including those ob- 
tained previously, (iv) Keep only noncquivalent graphs. 
The arrows in Fig. 2 indicate application of these rules to 
the second-order parent graph in the upper left corner. 14 
The same rules apply to the perturbation series for the 
"mass operator" V, except that in (i) the parent graphs 
are the diagrams for the total response function \. As 
discussed above, the introduction of V (or f xc ) is justified 
only if this function is free of singularities related to the 




FIG. 2: Examples of the second-order graphs for f xc . The ar- 
rows indicate an application of the rules (ii)-(iv) to the parent 
graph displayed in the upper left corner. 



FIG. 3: (a) General graph of the L-th order for the response 
functon x; (b) the diagram for the effective interaction which 
is produced from the general parent graph without self-energy 
insertions. 



KS e-h pairs in every order of the perturbation theory 
The analogy with the self-energy X does not, by itself, 
ensure this property, and we apply our graphical method 
to prove that this is indeed the case. 

Since the rules (ii)-(iv) deal only with the two-particle 
reducible graphs, a partial summation of the diagrams 
with the help of one- and two-particle irreducible ele- 
ments (self-energies and vertices) is possible. An exam- 
ple of a summation of all parent graphs with I vertex 
insertions (which divide each graph in I + 1 blocks) and 
with nik (k — 1, + 1) self-energies in every block is 
depicted in Fig. 3a. It is important that the diagrams 
generated by cutting two fermionic lines with the same 
frequency (e.g. in Fig. 1 and Fig. 2) have a similar struc- 
ture. They all describe scattering of KS particles by the 
xc potential. The sum of these graphs and the parent 
graphs is again a diagram of the same type as in Fig. 3a 
with the self-energy Eg = Gg 1 — G _1 . This definition 
of Ss ensures the equivalence of the KS and exact den- 
sities. The vertex T is defined in the usual way as the 
sum of the four-point functions which are irreducible in 
the e-h channel. After summation, the graph in Fig. 3a 
can be considered as a new parent graph which generates 
further diagrams for V via cutting only parallel e-h lines. 

Let us now consider the general parent graph of Fig. 3a 
at KS frequency = Ei — Ej, where Ej are the KS 
one-particle energies. This graph represents a L-th order 
correction (i.e. containing L irreducible elements) to the 
response function \- When the frequency u> approaches 
u>ij, the graph in Fig. 3a diverges. Integrating over inter- 
mediate frequencies in every internal block we find that 
the most divergent term behaves as l/(uj—uJij)( L+1 \ Ap- 
plication of our rules to the graph Fig. 3a generates a full 
set of diagrams for V{ui) i.e. the initial graph with only 
two external wiggled lines and the diagrams obtained by 
all possible cuttings of parallel e-h lines. All diagrams 
in this set have poles of various orders. However, they 
can be grouped in such a way that all singularities can- 
cel. The general proof is straightforward but lengthy and 
will be published elsewhere. Here we show this cancella- 
tion for a graph with L vertices, but with no self-energy 
insertions. 

Application of the rules (i)-(iv) to this graph leads 
to replacement of all internal KS two-particle functions 
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K s (uj,e) = G s (u J + e)Gs(s) by 

J(u,e,s') = K s (uj 7 s)8 e . £ , - K s {u},e)xs 1 {^)K s (uj 7 e'), 

where w and e are transferred (external) and internal 
frequencies respectively. The resulting contribution to V 
is shown graphically in Fig. 3b. The singularities at u> = 
ujij can potentially occur in the "resonant" part in every 
internal block Jjj(w, e, e'), which contains an electron in 
the state i and a hole in the state j. A summation over 
internal frequency gives for the "dangerous" contribution 

> Jy(o;,£,£ ) = xs M > 

^— ^ w — Uij oj — ujij u — Uij 

(4) 

where \ij) = ipi(r)i(j* (r 1 ) is the wave function of the res- 
onant e-h pair, \ij} — ipi(r)ipj (r) is the same function, 
but with equal coordinates of the electron and the hole 
and tpi(r) are the KS orbitals. Both terms in Eq. (4) are 
apparently singular. To show that these divergences can- 
cel, we single out the divergence in the KS susceptibility 
Xs = \ij}{ij\/(u-Uij)+Xr, where \r is the regular part. 
For the inverse Xg 1 (w) we have 

X S (U) ~ Xr {U) — — - (&) 

UJ-UJij +{lj\Xr 

Substitution of Eq. (5) to Eq. (4) gives a singularity-free 
result 



{ij\Xr 



(6) 



Similarly, zeroes of the external wiggled lines cancel the 
poles of the end blocks in Fig. 3b. 15 Thus the graph 
Fig. 3b is regular at KS resonances in every order of 
the perturbation theory, and V(u>) can be viewed as a 
mass operator similar to the one-particle self energy E. 
There is, however, one important difference. Whereas the 
self energy exactly reduces to the one-particle irreducible 
elements, the effective interaction V(co), even after can- 
cellation of singularities, still contains parts of the bare 
two-particle propagator with the resonant denominator 
being replaced by {ij\x r ~ 1 (u)ij)\ij} (see Eq. (6)). Since 
Xr 1 ( r : r ') m general goes to zero at |r — r'| — > oo and the 
functions \ij} have a normalization factor ~ 1/V, the 
matrix element {ij'IXr" 1 !^'} m ay vanish with the increase 
of the system size. For example, in a 3D semiconductor 
the inverse of the KS response function at ui = E g has an 
asymptotic behavior Xr 1 (Egi r ; r ') ~ l/l r — r 'l ( m this 
case Xr means the principal value of xs)- Therefore the 
matrix element \cv\Xr 1 (E g )\cv} vanishes as V~ x ^. This 
makes it problematic to construct approximations for f xc 
using perturbation theory, as is routinely done for E. 

The simplest of the perturbative approximations cor- 
responds to the first order in irreducible elements Eg 
and T (Fig. 4a), = X^X^Xs 1 , where X (1) is 

the first-order polarization loop. This class of approx- 
imations covers, for instance, the dynamic x-only OEP 
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FIG. 4: (a) Graphical representation of the general first-order 
approximation to the effective interaction; (b) self-energy 
Ss and irreducible vertex F for OEP; (c) Eg and T for RA ap- 
proximation. The double dashed line stands for the screened 
Coulomb interaction. 



(Ref. 9) and the Richardson-Ashcroft approximation 
(RA). 11 Voep(u) is given by the graph of Fig. 4a with 
Es and T taken in the first order in the Coulomb inter- 
action Vc (Ref. 8) in Fig. 4b. Similarly Vra(w) has the 
same form of Fig. 4a, but with the self-energy Eg and the 
vertex T as shown in Fig. 4c. OEP and RA are state-of- 
the-art approximations which accurately reproduce the 
correlation energy and plasma excitations of a homoge- 
neous electron gas. 12 However, both approximations may 
give uncontrollable results for the e-h excitation energies. 
In the following we show this by an explicit calculation of 
the second-order correction to the KS excitation energy 

U>ij. 



The first-order (with respect to Eg and T) correction 



Aw, 



(i) 



u $ ~ is s iven b y 



Alu 



(i) 



{ijiv^Mm^iijiwiij), 



(7) 



where the operator W(ri, r[; r 2 , r 2 ) is defined as follows 

W = EKn.ra^ri-i^J-E^rl.^n-ra) 

+ r efc (n,r , 1 ;r 2) r' 2 ). (8) 

The upper indexes in Eq. (8) indicate that the self- 
energies and the vertex are taken on the e-h mass shell. 
The equations (7)- (8) give a generalization of the x- 
only result, 8 which is equivalent to the first-order of the 
G6rling-Levy perturbation theory. 13 

Similarly we obtain the second-order correction 

+ fa1^ (1 H^j)|fcO("fc-"0{fc^ (1) K-)lu} ( p) 



where nk is the occupation number of the fc-th KS state 
andV^H is the effective interaction of the second order 
in Eg and T. If we ignore and iteratively solve 

Eq. (1) using only the first-order function (e.g. in 
OEP or RA approximation) in place of the kernel, the 
second-order correction Eq. (9) would contain only the 
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last two terms. The last term can be written as 



{ij\V( 1 \u ij ) X r(Ui j )VW(uJ ij )\ij} = 



_ cm; 



Uh2 



{ij\Xr 1 {^ )\iiY 



where we used Eq. (7) and Eq. (5). A similar calculation 
of the second term in Eq. (9) shows that it contains the 
same factor {ij\Xr Hence both terms in Eq. (9) 

which originate from and which would be the only 

contributions in OEP and RA, have exactly the same 
denominator as the higher-order diagrams after the can- 
cellation of the e-h singularities (see Eq. (6)). Clearly 
this denominator must appear in as well. The cal- 
culation shows that the first term in Eq. (9) indeed takes 
the form 

{y|V< 2 >K-)|y} = Ag } - 

- {ij\V( 1 \u ij ) X r(u ij )VW(u; ij )\ij}, (10) 



where 



,(2) 



E 



Wij - ion oiv 

(11) 

The second and the third terms in Eq. (10) cancel the sec- 
ond and the third terms in Eq. (9) and the second-order 



energy shift is simply given by Eq. (11). This also follows 
from a direct perturbativc solution of the BS equation. 
The cancellation found above is not accidental and can 
be proven in every order of the perturbation theory. 

Thus any finite-order approximation for V requires a 
consistent perturbative treatment of Eq. (1) to obtain 
a meaningful energy shift. In particular, the first-order 
approximations like OEP or RA are appropriate for the 
calculation of the excitation energies only in the first or- 
der in £5 and T. The naive second-order correction, 
(last two terms in Eq. (9)) and higher-order corrections 
are wrong and can even diverge in extended systems. A 
summation of all orders (exact solution of Eq. (1)) may 
produce a finite result, but it will contain uncontrollable 
errors. In recent TDDFT calculations for atoms it has 
been indeed observed that exact (in contrast to pertur- 
bative) solution of Eq. (1) does not necessarily improve 
the results. 5 
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matic rules for the xc potential v xc (Ref. 8). The only differ- 
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by far nontrivial fact leads to the conjecture that these 
rules should hold for any functional derivative 8 m E xc /8n m 
of the xc energy E xc . 

To simplify the derivation we assumed that the resonant 
transition u)ij is not degenerate. In the case of M-fold de- 
generacy, the denominator in Eq. (6) is replaced by the 
inverse of the M x M matrix {ij\Xr 1 l*'j'}> where \ij) and 
\i'j'} belong to the set of degenerate states. 



